Atmospheric Temperature vs Niño#
https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/doc/GHCND_documentation.pdf
import warnings
warnings.filterwarnings("ignore")
import os
import os.path as op
import sys
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import df2img
sys.path.append("../../../../indicators_setup")
from ind_setup.colors import get_df_col, plotting_style
from ind_setup.tables import plot_df_table
from ind_setup.plotting_int import plot_oni_index_th
from ind_setup.plotting import plot_bar_probs_ONI
plotting_style()
from ind_setup.core import fontsize
sys.path.append("../../../functions")
from data_downloaders import GHCN, download_oni_index
Define location and variables of interest#
country = 'Palau'
vars_interest = ['tmax', 'TMAX']
Get Data#
update_data = False
path_data = "../../../data"
Using Koror Station#
Analysis of how much the maximum and minimum temperatures over time are changing.
The analysis of the difference between these 2 variables will allow us to know how the daily variability is being modified
if update_data:
df_country = GHCN.get_country_code(country)
print(f'The GHCN code for {country} is {df_country["Code"].values[0]}')
df_stations = GHCN.download_stations_info()
df_country_stations = df_stations[df_stations['ID'].str.startswith(df_country.Code.values[0])]
print(f'There are {df_country_stations.shape[0]} stations in {country}')
if update_data:
GHCND_dir = 'https://www.ncei.noaa.gov/data/global-historical-climatology-network-daily/access/'
id = 'PSW00040309' # Koror Station
dict_min = GHCN.extract_dict_data_var(GHCND_dir, 'TMIN', df_country_stations.loc[df_country_stations['ID'] == id])[0][0]
dict_max = GHCN.extract_dict_data_var(GHCND_dir, 'TMAX', df_country_stations.loc[df_country_stations['ID'] == id])[0][0]
st_data = pd.concat([dict_min['data'], (dict_max['data'])], axis=1).dropna()
st_data['TMIN'] = np.where(st_data['TMIN'] >50, np.nan, st_data['TMIN'])
st_data['diff'] = st_data['TMAX'] - st_data['TMIN']
st_data['TMEAN'] = (st_data['TMAX'] + st_data['TMIN'])/2
st_data.to_pickle(op.join(path_data, 'GHCN_surface_temperature.pkl'))
else:
st_data = pd.read_pickle(op.join(path_data, 'GHCN_surface_temperature.pkl'))
ONI index#
https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php
p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
df1 = download_oni_index(p_data)
lims = [-.5, .5]
plot_oni_index_th(df1, lims = lims)
st_data_monthly = st_data.resample('M').mean()
st_data_monthly.index = pd.DatetimeIndex(st_data_monthly.index).to_period('M').to_timestamp() + pd.offsets.MonthBegin(1)
st_data_monthly
| TMIN | TMAX | diff | TMEAN | |
|---|---|---|---|---|
| DATE | ||||
| 1951-08-01 | 23.967742 | 31.045161 | 7.077419 | 27.506452 |
| 1951-09-01 | 23.993548 | 30.345161 | 6.351613 | 27.169355 |
| 1951-10-01 | 23.990000 | 30.673333 | 6.683333 | 27.331667 |
| 1951-11-01 | 24.509677 | 31.425806 | 6.916129 | 27.967742 |
| 1951-12-01 | 24.360000 | 31.023333 | 6.663333 | 27.691667 |
| ... | ... | ... | ... | ... |
| 2024-10-01 | 25.846667 | 29.926667 | 4.080000 | 27.886667 |
| 2024-11-01 | 25.855172 | 30.382759 | 4.527586 | 28.118966 |
| 2024-12-01 | 25.859259 | 30.148148 | 4.288889 | 28.003704 |
| 2025-01-01 | 25.788889 | 30.422222 | 4.633333 | 28.105556 |
| 2025-02-01 | 25.641176 | 30.170588 | 4.529412 | 27.905882 |
883 rows × 4 columns
rolling_mean = 6 #months
df1['tmin'] = st_data_monthly['TMIN'].rolling(window=rolling_mean).mean()
df1['tmax'] = st_data_monthly['TMAX'].rolling(window=rolling_mean).mean()
df1['tdiff'] = df1['tmax'] - df1['tmin']
df1['tmean'] = (df1['tmax'] + df1['tmin'])/2
fig, ax = plt.subplots(figsize = [15, 6])
df1.ONI.plot(ax = ax, color = get_df_col()[0], lw = 2)
ax2 = ax.twinx()
df1.tmin.plot(ax = ax2, color = get_df_col()[1], lw = 2)
# df1.tmax.plot(ax = ax2, color = get_df_col()[1], lw = 2)
# df1.tdiff.plot(ax = ax2, color = get_df_col()[1], lw = 2)
# df1.tmean.plot(ax = ax2, color = get_df_col()[1], lw = 2)
<Axes: >
low_lim = np.nanmin(df1.tmin)
fig, ax = plt.subplots(figsize = [15, 6])
df1.tmin.plot(ax = ax, color = get_df_col()[1], lw = 2)
ax.fill_between(df1.index, low_lim, df1.tmin, where = (df1.ONI > lims[1]), color = 'lightcoral',
alpha = 0.7, label = f'ONI over th: {lims[1]}')
ax.fill_between(df1.index, low_lim, df1.tmin, where = (df1.ONI < lims[0]), color = 'lightblue',
alpha = 0.7, label = f'ONI below th: {lims[0]}')
ax.fill_between(df1.index, low_lim, df1.tmin, where = ((df1.ONI > lims[0]) & (df1.ONI < lims[1])),
color = 'lightgrey', alpha = 0.075)
ax.legend(fontsize=fontsize)
ax.set_title('TMIN and ONI', fontsize = fontsize)
ax.set_ylabel('TMIN [°C]', fontsize = fontsize)
ax.set_xlabel('Time', fontsize = fontsize)
Text(0.5, 0, 'Time')
low_lim = np.nanmin(df1.tmax)
fig, ax = plt.subplots(figsize = [15, 6])
df1.tmax.plot(ax = ax, color = get_df_col()[1], lw = 2)
ax.fill_between(df1.index, low_lim, df1.tmax, where = (df1.ONI > lims[1]), color = 'lightcoral',
alpha = 0.7, label = f'ONI over th: {lims[1]}')
ax.fill_between(df1.index, low_lim, df1.tmax, where = (df1.ONI < lims[0]), color = 'lightblue',
alpha = 0.7, label = f'ONI below th: {lims[0]}')
ax.fill_between(df1.index, low_lim, df1.tmax, where = ((df1.ONI > lims[0]) & (df1.ONI < lims[1])),
color = 'lightgrey', alpha = 0.075)
ax.legend(fontsize = fontsize)
ax.set_title('TMAX and ONI', fontsize = fontsize)
ax.set_ylabel('TMAX [°C]', fontsize = fontsize)
ax.set_xlabel('Time', fontsize = fontsize)
Text(0.5, 0, 'Time')
plt.figure(figsize=(5, 4))
sns.heatmap(df1.corr(), annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('Correlation Heatmap')
plt.show()
df1.groupby(df1.index.year).max().tmean.min()
np.float64(26.77326692356587)
df1['tmean_ref'] = df1['tmean'] - df1.loc['1961':'1990'].tmean.mean()
df1['tmean_ref_min'] = df1['tmean'] - df1.groupby(df1.index.year).max().tmean.min()
df2 = df1.resample('Y').mean()
df2['oni_cat'] = np.where(df2.ONI > lims[1], 1, np.where(df2.ONI < lims[0], -1, 0))
fig = plot_bar_probs_ONI(df2, var='tmean_ref')
fig.suptitle('Temperature Anomaly over the 1961-1990 mean', fontsize = fontsize)
plt.show()
df_format = np.round(df1.describe(), 2)
fig = plot_df_table(df_format)
df2img.save_dataframe(fig=fig, filename="getting_started.png")